Black Hole Mergers and Unstable Circular Orbits 
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We describe recent numerical simulations of the merger of a class of equal mass, non-spinning, 
eccentric binary black hole systems in general relativity. We show that with appropriate fine- 
tuning of the initial conditions to a region of parameter space we denote the threshold of immediate 
merger, the binary enters a phase of close interaction in a near-circular orbit, stays there for an 
amount of time proportional to logarithmic distance from the threshold in parameter space, then 
either separates or merges to form a single Kerr black hole. To gain a better understanding of this 
phenomena we study an analogous problem in the evolution of equatorial geodesies about a central 
Kerr black hole. A similar threshold of capture exists for appropriate classes of initial conditions, 
and tuning to threshold the geodesies approach one of the unstable circular geodesies of the Kerr 
spacetime. Remarkably, with a natural mapping of the parameters of the geodesic to that of the 
equal mass system, the scaling exponent describing the whirl phase of each system turns out to be 
quite similar. Armed with this lone piece of evidence that an approximate correspondence might 
exist between near-threshold evolution of geodesies and generic binary mergers, we illustrate how 
this information can be used to estimate the cross section and energy emitted in the ultra relativistic 
black hole scattering problem. This could eventually be of use in providing estimates for the related 
problem of parton collisions at the Large Hadron Collider in extra dimension scenarios where black 
holes are produced. 



I. INTRODUCTION 

Gravitational physics is entering an exciting era. The 
construction of gravitational wave detectors that are ex- 
pected to be sensitive enough to observe many astro- 
physical phenomena where strong-field gravity plays an 
important role should teach us much about the cosmos 
and the structure of spacetime. Suggestions that we 
may live in a universe with more than 3 spatial dimen- 
sions 0, [|J offer the intriguing possibility that the Planck 
scale could be within reach of energies attainable by the 
Large Hadron Collider (LHC) [111, %. This implies 
that the LHC may be able to probe the quantum grav- 
ity regime, and that black holes could be produced in 
substantial quantities by the particle collisions. Simi- 
larly, cosmic ray collisions with the earth would produce 
black holes @, and this may be detected with current 
or near- future cosmic ray experiments 0]. Understand- 
ing the nature of black hole collisions within the theory 
of general relativity will be important in describing and 
interpreting many of these fascinating phenomena, if de- 
tected. 

The past couple of years has witnessed several break- 
throughs in numerical relativity 0, [f| [13], allowing for 
the solution of the field equations describing black hole 
mergers in many situations. However a vast region of 
interesting parameter space is unexplored, and it will be 
several years at least before a decent understanding of 
black hole collisions is achieved. In this paper we de- 
scribe simulations of the merger of a class of equal mass 
black hole binaries on initially eccentric orbits. The or- 
bits can be labeled with a one parameter family k loosely 
related to the initial tangential velocities of each black 
hole. We find intriguing behavior tuning k to what we 
call the threshold of immediate merger separating evolu- 



tions where the black holes either merge or do not dur- 
ing their first close encounter. The resultant evolution 
becomes exponentially sensitive to the initial parameter, 
and the binaries exhibit a period of "whirl" type behavior 
similar to that seen in geodesic motion [ll|, [l2| , orbiting 
rapidly in a near-circular configuration. Remarkably, sig- 
nificant amounts of gravitational radiation (w 1.0 — 1.5% 
of the rest mass energy per orbit) are still being emitted 
in this regime. Furthermore, based on the coordinate 
separation 1 , the binaries are orbiting within what would 
be the innermost stable circular geodesic of a Kerr space- 
time with angular momentum equal to that of the black 
hole that forms in the merger case. 

We show that the above behavior for equal mass 
binaries is analogous to evolution of similar classes 
of geodesies in a black hole background spacetime. 
Namely, if we define a one parameter family of equa- 
torial geodesies and tune to the threshold of capture, at 
threshold the geodesic will approach one of the unstable 
circular geodesies of the background spacetime, regard- 
less of whether the initial orbit is classified as elliptic or 
hyberbolic. Furthermore, if we calculate the instability 
(or Lyapunov) exponent of the orbits near the unstable 
circular orbits, we find numbers that are similar to that 
observed in this fully non-linear equal mass case. 

The binary black hole merger simulations described 
here are in the rest-mass dominated regime of the prob- 
lem. The close analogy with geodesic motion allows us 



1 Though this is not a gauge invariant quantity, though compar- 
isons between the extracted waveform and that estimated us- 
ing the quadrupole formula suggests the coordinate motion of 
the apparent horizons is quite well adapted to describing the 
situation UM- 
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to speculate about what might happen in the kinetic en- 
ergy dominated regime. This is relevant to the ultra- 
relativistic black hole scattering problem and thus might 
have application to the LHC if the Planck scale is below 
the maximum energies probed by the parton collisions. In 
particular, by finding the critical impact parameter and 
stability exponent of geodesic motion, estimating the en- 
ergy and angular momentum loss while the geodesic is in 
the whirl phase, and providing an estimate of the energy 
emitted for the head-on collision case in the full problem, 
one can come up with an estimate of the energy emitted 
to gravitational waves as a function of impact parameter. 
We also speculate that at threshold, all of the kinetic en- 
ergy of the system is converted to gravitational waves, 
which can be an arbitrarily large fraction of the total 
energy. 

The outline of the rest of the paper is as follows. In 
Sec|TT] we summarize the numerical code; in Sec lIIII we 
describe the simulation results; in Sec lIVI we describe the 
geodesic analog; in Sec|V]we discuss all the results, spec- 
ulating about what may happen beyond the equal mass, 
non-spinning case, and how the results might carry over 
to the kinetic energy dominated regime and be applied 
to the BH scattering problem; in Sec|Vl]we provide some 
concluding remarks and a discussion of future related 
work. The technique we use for geodesic integration is 
described in the appendix. 
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In the above, T c ab are the Christoffel symbols 
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T a b is the stress energy tensor with trace T, H a are the 
GH source functions defined via 



(4) 



k is a parameter multiplying the constraint damping 
terms, n a a timelike vector, and C a are the constraints: 



C a = H a - Dx a . 



(5) 



Any solution to the Einstein equations must have C a = 0; 
in a consistent numerical evolution the constraints will 
be zero to within the truncation error of the numerical 
scheme. For n a we use the usual hypersurface unit nor- 
mal vector n a = l/a(d/dt) a , where a is called the lapse 
function. 

We couple in a massless scalar field $ as the matter 
source, which satisfies the wave equation 



II. OVERVIEW OF THE EQUATIONS AND 
SOLUTION METHOD 



The Einstein Field Equations (EFE) in generalized 
harmonic form have been discussed in much detail else- 
where [3 [TElilM EES SI Si 11, [M SIM S3, 

SIJ so for brevity here we simply list the equations and 
briefly mention the numerical code solving these equa- 
tions. 



A. formalism 

We solve for a spacetime described by the line ele- 
ment ds with metric tensor g ab and coordinates x a = 

(t,x,y,z) 2 



□$ = 0, 
and has a stress energy tensor 

T ab = 2$ $ b - g ab $, c $' 



(6) 



(7) 



The following equations are used to evolve the source 
functions: 



a - 1 

a" 



&H ttU n u , Hi = 0. (8) 



ds 2 = g a bdx a dx , 



(1) 



This equation for Ht is a damped wave equation with 
a forcing function designed to prevent the lapse a from 
deviating too far from its Minkowski value of 1. The pa- 
rameter £2 controls the damping term, and £1 , 77 regulate 
the forcing term. The particular parameter values are as 
described in 12111. 



B. Boosted Scalar Field Initial Data 



using the EFE in generalized harmonic (GH) form with 
constraint damping [29l . |3C| : 

-^g cd g a b,cd + g cd (, a gb)d, c + H( a ,b) ~ HdT d ab 



2 We use a comma (,) to denote partial differentiation, and we use 
units where Newton's constant G = 1 and the speed of like c = 1. 



For initial conditions we use two boosted scalar field 
pulses, as described in detail in [2l{ . These pulses very 
quickly undergo gravitational collapse and form a pair of 
black holes in a bound orbit, with most (« 85%) of the 
scalar field energy falling into the black holes, the rest ra- 
diating away on roughly the light-crossing time scale of 
the orbit. At the initial time we assume the correspond- 
ing spatial geometry is conformally flat and maximal, and 
solve the Hamiltonian and momentum constraint equa- 
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tions together with the maximal slicing condition for self- 
consistent initial conditions. For the remaining coordi- 
nate degrees of freedom we choose the spacetime slice to 
be harmonic at t = 0. 



C. Numerical code 

The evolution code, described in detail in 1!) . uses 
second order accurate finite difference discretization with 
adaptive mesh refinement, a coordinate system compact- 
ified to spatial infinity and excision techniques to deal 
with the singularities inside of the black holes. Proper- 
ties of the black holes are measured using apparent hori- 
zon (AH) properties, and gravitational wave information 
is extracted using the Newman Penrose formalism with 
the extraction radius placed a distance of 50m from the 
origin, where m = mi + m 2 is the total, initial AH mass 
of the individual black holes measured after the majority 
of scalar field energy has been accreted. 



III. THE THRESHOLD OF "IMMEDIATE" 
MERGER 



Here we present one of the main results of this paper, 
namely we give evidence that there are regions in the 
parameter space of the two body problem in full gen- 
eral relativity where the black holes evolve toward an 
unstable circular orbit, remain in that configuration for 
an amount of time sensitively related to the initial condi- 
tions, then either plunge toward coalescence or separate. 
In the case where the black holes separate after the circu- 
lar motion they could possibly merge at some time in the 
future. These regions in parameter space can be found 
by examining suitable one parameter families p of initial 
conditions, where pi < p < p s smoothly interpolates be- 
tween an evolution with p = pi where a merger occurs 
promptly within t = U, and an evolution with p = p s 
where after some amount of time t s 3> £j no merger has 
occurred. The unstable circular orbit is approached near 
the threshold of immediate merger at p p*, where for 
p < p* merger occurs promptly, while for p > p* it does 
not 3 . The number of orbits n observed near threshold is 
found to scale as 



oc \p-p*[ 



(9) 



Note that due to the energy loss via gravitational radi- 
ation the threshold cannot be "sharp" , i.e. if the time 
tm(p) to merger is plotted as a function of p, this will 
not have a discontinuous step at p — p* . There will be a 



3 Note that depending on the parameter that is being varied the 
magnitude of p relative to p* denoting prompt merger can be 
inverted. 



maximum number of orbits N for a given class of initial 
conditions, and what from a distance might appear like 
a step function will be resolved into a smooth transition 
over a region of size Sp rs e~ N / J . 

We cannot prove some of the statements made in the 
preceding paragraph, in particular since the full numeri- 
cal simulations are so computationally expensive we have 
only studied one class of initial conditions in detail. How- 
ever, these simulation results, the striking similarity be- 
tween them and the geodesic analogue presented in the 
next section, and trying to understand what must hap- 
pen at a generic threshold as discussed in Sec fVl provides 
quite a compelling case for this behavior. 



A. Simulation Results 

Here we describe results from the evolution of a class 
of scalar field collapse binary (SFCB) simulations dis- 
cussed before in 21| . The new simulation data presented 
here includes several higher resolution runs tuned closer 
to the threshold of immediate merger, and so now we 
have significantly more confidence that we are converg- 
ing to this phenomenon in the two body problem. Below 
we give a brief summary of the initial conditions, and 
focus on evolution results of relevance to classifying and 
understanding the immcdiatc-merger-threshold scenario. 
More background details can be found in 21]. 

For initial conditions we begin with two identical 
boosted scalar field distributions, one located at a coordi- 
nate location of (x, y, z) — (4.45m, 0, 0) and given a boost 
with boost parameter k in the positive-?/ direction, while 
the other is located at (x, y, z) — (—4.45m, 0, 0) and given 
a boost k in the negative-y direction (the proper separa- 
tion is initially 10.8m). The scale m here is the sum of 
apparent horizon masses measured around t — 15m after 
evolution has begun, which is after essentially all of the 
collapsing scalar field energy has either accreted into the 
newly formed black holes or is on its way to escaping to 
infinity. Approximately 85% of the initial scalar field en- 
ergy falls into the black holes. Thus k labels our family 
of initial conditions. Note that k is related to though not 
exactly the same as the initial velocities the black holes 
will have. With k = we have a head-on collision, i.e. 
a prompt merger, while for k sufficiently large the black 
holes arc deflected but fly apart. Thus k describes an 
appropriate family to study the threshold of immediate 
merger. 

Fig. Q] shows n(k) © for the cases that merged (A; < 
k*), while Fig. [5] shows n(k) for the cases that separated 
after the initial whirling motion. In the former plot n 
is calculated as the total phase angle <f> divided by 2ir 
that one of the black holes moves through before a com- 
mon AH is detected, while in the later case is the total 
4> evolution divided by 2n undergone by the black hole 
in returning to its original coordinate distance from the 
origin. An example of the orbital trajectories is shown in 
Fig. O Fig. 2] shows the total energy radiated in grav- 
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FIG. 1: The number of orbits n versus logarithmic distance 
of the initial boost parameter k from the immediate merger 
threshold k* , for evolutions that did result in a merger. Re- 
sults from the three resolutions are plotted, and a least- 
squares fit to each set of data assuming the relation Q. 



FIG. 2: Data as in Fig. [TJ though here from evolutions that 
did not merge during the time of the simulation (i.e. k > k"). 



itational waves from these simulations, Fig. [5] shows a 
sample of the waveform from a couple of the simulations 
measured using both the Newman-Penrose scalar vE^ and 
the quadrupole formula, and Fig. [5] shows the energy flux 
(calculated using ^4) from four sample simulations. 

From the data shown in FigQ]we estimate 7 = 0.35 ± 
0.03 for this family of initial data. It is difficult to calcu- 
late the uncertainty in this quantity. In theory conver- 
gence testing should be sufficient, though here each res- 
olution has a different number of points, so the intrinsic 
error in a linear regression fit will be different. Further- 
more, we are assuming @ holds, and it most certainly 
does not exactly. Also, since each set of simulations span 
a different range in k — k* this will cause some variation 
in the measured 7 in addition to truncation and small 
sample size errors. We have therefore simply taken the 
uncertainty to be the difference in 7 between the medium 
and higher resolution simulations (which is roughly what 
it would be if only truncation errors were responsible for 
the differences). A summary of the three characteristic 
resolutions used is listed in Table HI 



IV. THE GEODESIC ANALOGUE 

To gain insight into what is happening when tuning to 
the threshold of immediate merger it is useful to compare 
to a geodesic analogue of this behavior. Specifically, we 
will play the same threshold game with geodesies in a 
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FIG. 3: Plots of the orbital motion from the two higher reso- 
lution simulations (h/2) tuned closest to threshold (only the 
coordinate motion of a single black hole — initially at positive 
x — is shown for clarity). The dashed curve is the case result- 
ing in a merger, and the curve ends once a common apparent 
horizon is first detected, while for the solid curve the black 
holes separate again and here the curve ends when the simu- 
lation was stopped. 
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FIG. 4: The total energy radiated in gravitational waves plot- 
ted as a function of logarithmic distance from the immediate 
merger threshold. We have overlayed the data from both su- 
per and sub critical cases, though for clarity have only added 
the linear regression line for the cases that merged. 



"Resolution" 


wave- zone res. 


orbital-zone res. 


BH res. 


h 


1.7m 


0.23m 


0.057m 


3/4 h 


1.3m 


0.17m 


0.043m 


1/2 h 


0.85m 


0.12m 


0.029m 



TABLE I: The three sets of characteristic resolutions used 
in simulation results presented here, where each resolution 
is labeled relative to the coarsest resolution h. The grid is 
adaptive with a total of 8 levels of refinement, and the coor- 
dinate system is compactified. The wave zone is defined to 
be at r = 50m, the orbital zone within about r = 10m and 
the black hole zone is within 2 — 3m of each apparent hori- 
zon, where m is the initial sum of apparent horizon masses, 
measured at t = 15m to account for early scalar field accre- 
tion. At this time the coordinate radius of each apparent 
horizon is w 0.52m; thus for the highest (lowest) resolution 
case roughly 36 (18) grid points span each horizon. A CFL 
(Courant-Friedrichs-Lewy) factor of 0.15 was used in all cases. 



Kerr background by constructing a one parameter fam- 
ily p of geodesies where for p < p* the geodesies even- 
tually fall into the black hole, while for p > p* they do 
not. For families of bound (elliptic) orbits, the latter sub- 
set of parameter space exhibit zoom-whirl behavior — the 
geodesies start some distance from the central black hole, 
move in close to the black hole where they whirl around 
for several orbits, then zoom out again to the original 
distance; this behavior then repeats. For the case p < p* 
the initial behavior is similar, namely the geodesic moves 
in from a distance and starts whirling about the black 
hole, but then plunges into the black hole. The num- 
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FIG. 5: The gravitational waves emitted during a merger 
event. Here we show the real part of the dominant spin weight 
-2, £ = 2, m = 2 spherical harmonic component of ^4. The 
solid curve is the merger case tuned closest to threshold (from 
the higher resolution simulations) . The dotted curve was com- 
puted by taking the coordinate motion of the AH's from this 
simulation (see Fig(3J, assuming they represent the locations 
of point particles of mass m/2, and plugging the result into 
the quadrupole formula; this waveform ends once a common 
horizon is formed. The dashed curve shows the same infor- 
mation as the dotted curve, but here the data is from the 
non-merger case tuned closest to threshold. Note that the 
two curves from the quadrupole formula were shifted in time 
to account for the propagation time to the sphere at r = 50m 
where ^4 was measured. 



ber of whirl orbits increases as p approaches p* , going to 
infinity in the limit. In this limit, the whirl part of the or- 
bit comes arbitrarily close to one of the unstable circular 
geodesies of the spacetime (we are restricting attention 
to equatorial orbits here). Qualitatively the same behav- 
ior is seen in families of unbound (hyperbolic) geodesies, 
the only difference is that there is only one episode of 
whirling before the subset of geodesies that do not fall 
into the black hole escape to infinity. 

This behavior is not only qualitatively similar to what 
is observed in the full merger simulations, there is also 
quantitative similarity comparing the scaling exponent 
7 in ([9]). Note that there even is a geodesic analogue 
is surprising, in particular if we take what the analogue 
is telling us at face value: in this fine tuned regime the 
binaries are approaching an unstable circular orbit. In 
the full problem, one might guess that even if the bina- 
ries do temporarily approach what is an unstable circu- 
lar orbit, then surely the "perturbation" implied by the 
rather strong gravitational wave emission with this close 
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FIG. 6: The energy flux in gravitational waves from four (h/2 
resolution) simulations; three resulted in mergers, the fourth 
(long-dashed curve) not. 



separation would quickly force a merger? However this 
does not seem to happen, and moreover we speculate that 
with sufficiently fine-tuned initial conditions the binary 
can remain in this regime until close to as much energy as 
is theoretically possible is radiated away in gravitational 
waves (more on this in Sec fV} . 

In the remainder of this section we describe the 
geodesic analogue in detail, first giving an analytical cal- 
culation of the scaling exponent 7 in Sec lIV Al then show- 
ing results from numerical integration in Sec lIV Bl 



A. Calculating 7 from perturbation theory 

Here we calculate 7 by finding the instability exponent 
A of unstable circular geodesies. Our analysis mirrors the 
technique and formalism of 31] for Schwarzschild space- 
time; here we extend the result to equatorial orbits of 
Kerr. We will perform the calculation in Boyer-Lindquist 
coordinates, in which the Kerr metric takes the following 
form 



ds 2 = 



2mr\ , 9 S , 9 

r + A dr 



on n 4mar sin 2 9 
-R 2 sin 2 9d4> 2 (tydt, 



where 



r 2 + a 2 cos 2 1 



(10) 



(11) 



R 



2ma 2 r sin 2 9 



A = r 2 + a 2 — 2mr. 



(12) 
(13) 
1 the 



m is the total mass of the black hole, and J = 
total angular momentum. 

First, we want to relate A, the growth rate of radial 
perturbations of the orbits in coordinate time t, to 7, that 
characterizes how the number of orbits increases as the 
logarithmic distance to the critical parameter decreases. 
The orbital angular frequency is 



U = (f>, 



(14) 



where the overdot (') denotes d/dt. We will start with a 
geodesic with initial conditions r(t = 0) = r^ + Sr^, where 
ro is the radial coordinate of one of the unstable circular 
orbits, and Sro is a small perturbation. We denote the 
corresponding constant orbital frequency of the circular 
orbit with u>o, which is 



UJo 



ma ± \fi 



(15) 



where + sign is for co-rotating, the — sign for counter- 
rotating geodesies. To linear order the growth of the 
perturbation Sr(t) = r — ro is given by 



5r(t) = 5roe 



m 



(16) 



Taking the absolute value and then natural logarithm of 
both sides we get 



ln\Sr(t)\ =ln\Sr \ + Xt. 



(17) 



Perturbation theory breaks down when 5r(t) ~ 1, which 
is also the time when the geodesic will either leave the 
vicinity of the black hole, or fall into it. By this time the 
number of orbits that have been completed is n — ujQt /2ir. 
Finally, we note that Sro will be proportional to p — p* 
for a family of geodesies that approach this unstable orbit 
when p = p*. Substituting all this into (jTTJ) and solving 
for n gives 



In \p-p*\-^\ + const., 



from which we can read off 7 ©: 



7 



^0 
2ttA' 



(18) 



(19) 



Now we turn to the calculation of 7. We restrict at- 
tention to equatorial (9 = 7r/2) geodesies, for which the 
corresponding Lagrangian is 



2m\ 'o 1* 2 t n 4ma > 
2£ = -(l ] t 2 + —r 2 + R 2j - 2 * 



cf>t, (20) 



where R 2 , 



2 (1 + 2m /r), and the prime (') de- 



notes differentiation with respect to affine parameter s. 
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The momentum p q conjugate to geodesic coordinate q 
is SC/5q , and in terms of p q and q the Euler-Lagrange 
equations for each pair (p q ,q) is 



dp 
ds 



5q 



(21) 



Given that £ does not explicitly depend on t and <f> we 
immediately obtain the following two first integrals of 
motion, which as usual we identify as the energy E and 
angular momentum L of the geodesies: 



E = - Pt = 
L = p^ = Rl 



2m\ > 2ma i , . 

1 t + , 22 

r I r 



Iraa 



-t . 



(23) 



A third constant of motion comes from the normalization 
condition for timelike geodesies 



2m\ / 9 r > „ 9 4ma ,/ / . „. 

1 )t +-r r 2 + Rl4> 2 4>t =-1. 24 

r I A r 



In principle we can use the above three equations to find 
a first order differential equation for r(s), though instead 
we will follow [3l| and evolve the pair (r,p r ) using (|2"Tj) . 
We shall also explore the dynamics in terms of coordinate 
time t rather than affine parameter s; the transformation 
between them can be obtained from (|22I23| : 



t 



1 

A 



(ERl ~ 2maL/r 



The resultant system of equations is 

A p r 

r = 

rH' 

t 



p r = — [u> 2 (r 3 — ma 2 ) + m(2au> — 1)] 

Pr 

t'r 3 



Pr r 2 1 

[a — mrj 



Following [3l| , we write the above as 



dXjjt) 
dt 



= Hi(Xj) 



(25) 



(26) 



(27) 



(28) 



where Xi(t) — (r(t),p r (t)) and Hi(Xj) are the right hand 
sides of the corresponding equations. We linearize the 
equations about circular orbits, namely let Xi(t) = Xio + 
SXi(t) where X iQ = (r = r Ql p r = 0), and only keep terms 
linear in 5Xi(t): 



d6Xi(t) 
Jt 



KaSXAt), 



where 



(29) 



(30) 



The quantity A that we are interested in is the positive 
real eigenvalue of (f30|) . which exists for the range of ?'o 
corresponding to unstable circular orbits. We can then 
substitute this into (|19[) to find 7. After a tedious but 
straight-forward calculation we get 



„2 r 



7 



2tt 



4to 



3r 2 A H {rRluo 2 - Amauj - r + 2m) 



-1/2 



(31) 

where in the above u> and r are evaluated at luq and tq 
respectively. Note that although the preceding calcula- 
tion is gauge dependent, the final result (f3i"j) only directly 
refers to the radial coordinate r, and indirectly to the az- 
imuthal coordinate <fi in that we measure the number of 
orbits n completed by the geodesic, r can be eliminated 
entirely from the expression in favor of the proper cir- 
cumference 27ri?o of the geodesic (though for simplicity 
we have not done so), and given the axisymmetry of the 
spacetime one can unambiguously define n. 



B. Equatorial geodesic orbits in Kerr 

In the preceding section we calculated 7 by examin- 
ing the growth of a perturbation of an unstable circular 
geodesic. However, it might not be obvious that one pa- 
rameter families of geodesies that smoothly interpolate 
between capture and non-capture will approach one of 
these unstable orbits as the limiting case between cap- 
ture and non-capture. Here we show a few examples that 
this is indeed the generic behavior at threshold for equa- 
torial geodesic families, though the particular unstable 
circular orbit that is approached depends on the initial 
conditions. 

Figs. [7] and [5] show a couple of examples for two fine- 
tuned families of initial data in Schwarzschild spacetime 
(a = 0), integrated using the method explained in Ap- 
pendixlAlfand note that we have integrated the geodesies 
in Cartesian coordinates in the Kerr-Schild form of the 
metric, which is horizon penetrating). Fig. [7] is an ex- 
ample of a hyperbolic family of orbits, Fig. [5] an elliptic 
family. This again illustrates that no special care need be 
taken in choosing a class of orbits to exhibit this behav- 
ior, merely that the orbits can be labeled by a parameter 
that smoothly intersects the threshold of capture. Fig. [5] 
shows 7 as a function of ro, the radius of the geodesic in 
the whirl-phase, for a range of values of the spin param- 
eter a. The lines in the figure were calculated using (|31[) . 
and for three representative cases of a we overlay the 
results from a numerical calculation based on the same 
method used to calculate 7 for the fully non-linear prob- 
lem described in Sec lIIIl Note that each point from the 
numerical calculation represents fine-tuning a family of 
geodesies to the threshold of capture, and we typically 
tuned the initial conditions to within 1 part in 10 12 (be- 
ing geodesic integration on a fixed background this is a 
"cheap" problem) . It was also not difficult to find sets of 
one parameter families that at threshold spanned most 
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of the range of unstable circular orbits. For example, 
repeating the exact numerical setup of the binary black 
hole problem, we can choose the geodesic to have some 
initial coordinate (x, y, z) — (d, 0, 0) and an initial veloc- 
ity (0,Vj,,0). Keeping d fixed and tuning v y to thresh- 
old we approach an unstable circular orbit the radius of 
which depends on d. So repeating this threshold search 
for a range of d we can map out most of the range of 
unstable circular orbits. 

On Fig. [9] we also plot an ellipse representing 7 and tq 
measured in the equal mass problem described in Sec lIIH 
the size of the ellipse depicts the numerical uncertainty 
in this "point" . For ro we have used the coordinate sepa- 
ration between the black holes, which, based on how well 
this quantity reproduces the gravitational waves emitted 
when plugged into the quadrupole formula (see FigfS]) 
suggests this is a reasonable distance indicator not sub- 
ject to severe gauge artifacts 4 . The region of the diagram 
where 7 and ro from the binary merger simulations fall 
lends support to the idea that the geodesic analogue is 
describing what is happening in the full problem — the co- 
rotating geodesic orbit with the same value of (^0,7) as 
the non-linear problem gives a spin parameter of ss 0.5, 
which is not too far from that of the final merged black 
hole of « 0.7. 



V. DISCUSSION 

In Sec lIIII we gave evidence of a "threshold of immedi- 
ate merger" in the binary black hole problem for one pa- 
rameter families of initial data that interpolate smoothly 
between a scenario of prompt merger at one extreme 
and non-merger (or merger much further into the future) 
at the other extreme. Near threshold evolutions were 
marked by a period of near-circular orbital evolution at 
very close separation, significant gravitational wave emis- 
sion, and exponential sensitivity to initial conditions. In 
Sec lIVI we demonstrated that this behavior was quanti- 
tatively similar to equatorial geodesic behavior in a Kerr 
background, where the corresponding threshold is that 
of capture of the geodesic. Without the geodesic ana- 
logue this behavior in the full problem might seem very 
bizarre, and indeed even with the geodesic analogy it is 
still rather surprising. For given how unstable the cir- 
cular orbits are that are approached at threshold, one 
would imagine that a "perturbation" such as the emission 
of a significant amount of gravitational radiation per or- 
bit would make such an orbit unattainable away from the 
geodesic approximation. However, as we will try to argue 
in Sec lVAl below, this behavior is almost "obvious" if we 
try to imagine the possible orbits that could arise in the 
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4 Another suggestion might be to use the proper distance between 
the two horizons; however this, when used in the quadrupole 
formula, overestimates the energy by almost a factor of 2. 



FIG. 7: Plots of two timelike geodesies of Schwarzschild 
(m — 1) spacetime tuned close to the threshold of capture. 
The geodesic integration was started at (x, y, z) = ( — 10, 0, 0), 
the geodesies given a velocity of v, and the parameter used 
to find the threshold was the initial angle 9 of the tangent 
vector relative to the x-axis. Here, the two geodesies evolve 
toward an unstable circular orbit at 3.54, remain close to it 
for roughly 6 orbits, then the one with 9 < 9* (dotted curve) 
falls into the black hole while the other (solid curve) escapes. 
The dashed circle is the location of the event horizon. Note 
that the unstable circular orbit that is approached has r < 4, 
which is consistent with the orbit being unbound (see for ex- 
ample [3^| )■ 



full problem. The argument applies to the generic sce- 
nario of unequal mass, arbitrary spin initial conditions, 
however the argument cannot claim that generically there 
will be any relationship with unstable geodesies of black 
hole spacetimes. Near-threshold orbits almost certainly 
do not have astrophysical significance due to the fine- 
tuning required to reach them, however if large extra 
dimensions exists and will result in significant black hole 
production at the LHC, this behavior could have some 
application there — we discuss this in Sec. IV Bl 



A. Beyond equal mass non-spinning mergers 

Consider a smooth one-parameter family b of initial 
conditions for two black holes of arbitrary masses and 
spins that straddle a regime of immediate merger. We 
further assume that the parameter is monotonic in that 
the threshold is only crossed once as b is varied (this also 
assumes that the threshold is not fractal in nature). Let 
b = b± denote one extreme (merger), b — 62 the other 
(deflection). Now choose a value of b halfway between 
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FIG. 8: A second example of timelike geodesies of 
Schwarzschild (m = 1) spacetime tuned close to the thresh- 
old of capture. Here, the geodesic integration was started at 
(x,y,z) — (—4.5,0,0), and the parameter tuned to thresh- 
old was the initial velocity v = (0, v y , 0) in the y direction. In 
other words, we are effectively starting on an unstable circular 
orbit with a "perturbation" v y — Vy. The two geodesic orbits 
shown here remain close to r = 4.5 for about 8 orbits, one 
(dotted curve) then falls into the black hole, while the other 
escapes. Since r > 4 this is a bound orbit (see for example 
[3^]). and the geodesic "zooms" out to some distance before 
returning to r m 4.5 to undergo another "whirl" episode. Here 
we show two full zoom-whirl episodes (and note that this is 
not a closed orbit; the integration was stopped during the 
start of the third whirl phase). 



these two: 63 = (b\ + &2V2. By the assumed smoothness 
and monotonicity of b the resultant orbit must lie between 
the first two orbits , and given that one orbit is anchored 
in a merged black hole, the black holes in the new orbit 
must spend more time in close proximity than either of 
the preceding cases before merging or separating — this is 
illustrated in Fig. [101 Continuing the bisection therefore 
forces the corresponding binary system to approach an 
every lengthier whirl-like period of evolution. However, 
since more time is spent in a tight orbit, more energy will 
be radiated, and if the binary was initially unbound, after 
a certain amount of fine-tuning the system will become 
bound. Continuing the process beyond this point (or 
from initial conditions that were bound to begin with), 



5 Though with one or both of the black holes possessing spin the 
orbital plane will in general not lie on a plane, making the tra- 
jectories curves in 3D space and preventing a simple and precise 
notion of "between". 




FIG. 9: Plots of 7(7*0) for Kerr equatorial geodesies (|3ip . for 
values of the spin parameter a ranging from a — —0.9 (right- 
most curve) to a = 0.9 (left-most curve) with intervals of 0.1. 
For three values of a (0, ±0.7) we plot 7(ro) calculated us- 
ing geodesic integration tuning to the threshold of immediate 
capture. Specifically, each point depicted in the figure shows 
results of a bisection search to the threshold of capture for a 
one parameter family of geodesies; the value of ro is extracted 
from the nearest-to-threshold solution, and 7 is obtained from 
data similar to that shown in Fig[T] The dashed ellipse shown 
in the figure is the value of 7 (enlarged from a point to include 
estimated uncertainties) extracted from the fully non-linear 
evolutions described in Sec lIIII This is in rather remarkable 
agreement with the geodesic results given that the final spin 
of the black hole is 0.68 ± 0.05, and the trend from the sim- 
ulations shows ro is slowly decreasing as one approaches the 
threshold, suggesting that the ellipse would shift to the left if 
we could tune closer to threshold. The near-critical geodesies 
at the center of the dashed ellipse have an eccentricity of 
e ~ 0.5, which might be one way to define what the effective 
initial eccentricity of the equal mass system is. 

the effective apoapsis of the orbit for the cases that sepa- 
rate must decrease as the threshold is approached. Grav- 
itational radiation therefore puts a limit to this process, 
ending it once enough energy has been radiated so that 
the next apoapsis decreases to the radius of the orbit the 
binary is on during the whirl phase. Of course, approach- 
ing this point the notion of an "immediate merger" will 
become ill-defined, for ever sooner after separating the 
binaries will merge and thus it will become impossible to 
differentiate between prompt merger and separation. 

The preceding argument appeals to continuity of the 
orbital trajectories as a function of a smooth initial data 
parameter b. One way of preserving continuity, yet hav- 
ing the whirl regime terminate before all possible energy 
has been radiated, is if at some distance from threshold 
the encompassing black hole forms while the two black 
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FIG. 10: A schematic illustration of why an immediate thresh- 
old must exist for the generic scattering problem. Each panel 
shows several hypothetical trajectories of one black hole in a 
two body interaction in the center of mass frame; for clarity 
the second black hole is not shown. The two thicker solid 
curves represent impact parameters bracketing the threshold, 
i.e. in one case there is a merger, in the other a deflection. The 
two thinner dashed curves represent the two possible trajec- 
tories for an impact parameter between the bracketing cases. 
Going from the top to bottom panel we are successively tuning 
closer to threshold. It is difficult to imagine how a bisecting 
trajectory could not follow a path between the two bracket- 
ing trajectories as illustrated, and so tuning to threshold one 
forces the binaries into a close "whirling" configuration. 



holes are still in the whirl phase. However, given that 
both the simulations shown in Sec lIIII and the geodesies 
in Sec|W]have a distinct plunge phase for the cases that 
do merge, we think this possibility is unlikely. One inter- 
esting question then is how much binding energy could be 
extracted at threshold via gravitational waves. If cosmic 
censorship holds then an upper limit is given by Hawk- 
ing's area theorem. The exact amount depends on the 
angular momentum of the final black hole; if it is non- 
rotating a maximum of 29% of the rest mass could be 
radiated; for the case studied here if the trend continues 
and the final spin remains a « 0.7, 24% could be ra- 
diated. In the high-speed scattering problem where the 
kinetic energy of the black holes dominate the net energy, 
if the argument outlined in the previous paragraph holds 
then it may be possible for all the kinetic energy of the 
system to be emitted as gravitational waves at thresh- 
old. This can be an arbitrarily large fraction of the total 
energy of the system. 

We note that for the above arguments to be valid it is 
crucial that two distinct end-states exist in the scattering 
problem. In Newtonian gravity for two point masses, 
for example, this type of behavior cannot exist because 
there is no merger end-state. It is also interesting that in 
the geodesic problem the whirl-type behavior that arises 
at threshold is intimately connected to the existence of 
unstable circular orbits; this makes it tempting to think 
that there may be some deep connection between the 
existence of unstable geodesies in Kerr and the fact that 
black holes merge when they collide in general relativity. 



B. Black hole scattering at the LHC 

One interesting application of this apparent analogy 
between the threshold of capture of geodesies and the 
threshold of immediate merger in the full problem is to 
obtain some understanding of the black hole scattering 
problem. This is of interest in the search for black holes 
that might be formed at the LHC (for a few review arti- 
cles see 0, [HI, H3| ) . Consider the scattering of two black 
holes colliding with an impact parameter b. We would 
like to know E r (b) 1 the energy radiated as a function of 
the impact parameter 6, and the threshold impact param- 
eter b* below which a merger results (this would corre- 
spond to the threshold of black hole production in a parti- 
cle collision of sufficiently high energy where the classical 
general relativistic description of the event holds). Be- 
low we illustrate how to find an approximation to E r (b) in 
the 4 dimensional case; certainly to be of relevance to the 
LHC this calculation will need to be carried out in higher 
dimensions, and many of the approximations made could 
be improved upon, though here we are more interested in 
presenting the main ideas. Note that higher dimensional 
black hole spacetimes in general do not posses stable cir- 
cular orbits [35l [3rj \3l\ , though what is crucial here is the 
existence of unstable circular orbits. 

For concreteness we consider the interaction of two 
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equal mass, non spinning black holes of total rest mass 
77i traveling with large center-of-mass speeds v. In this 
regime the total energy of the system E = Tm = 
m/yl — v 2 is dominated by the kinetic energy of the 
black holes. Also, for sufficiently large T any spins or elec- 
tric charges of the black holes will become irrelevant, as 
these quantities are Lorentz invariant and thus the ratio 
of spin-to- kinetic and/or electromagnetic-to- kinetic en- 
ergy of the system goes to zero. So the non-spinning, un- 
charged case will be a good approximation to the generic 
high energy scattering problem 6 . Denote the threshold 
of immediate merger by 6 = b*. Guided by the results 
presented in earlier sections, we will assume the following 
key pieces of information to obtain an approximation to 
E r (b): 

• evolutions near the threshold of immediate merger 
b* are characterized by the scaling relation ([9]): 

e™ oc \b- 6*| 7 

• the critical impact parameter b* and the scaling 
exponent 7 can be found by considering the analo- 
gous problem of geodesic scattering in a Kerr back- 
ground. The ADM mass of the background space 
time is E (which sets the scale for 6), and has a 
spin parameter a equal to that of the black hole 
that forms in the full problem for b slightly less 
than b* . 

• during the whirl phase, a constant fraction e of the 
remaining energy E — E r (t) of the system is radi- 
ated per orbit 

There are a couple of additional bits of information that 
will be needed to complete the calculation — the energy 
E r o = E r (0) emitted in a head-on collision, the value of 
a for the background Kerr spacetime, and the fraction 
e of energy radiated per orbit. We will use an existing 
estimate of E r o from the literature, argue below that a 
slightly less than extremal (a = 1) is the relevant pa- 
rameter in the large T limit, and use the quadrupole for- 
mula for e. Using the quadrupole formula together with 
geodesic motion is similar in spirit to the "semirelativistic 
approximation" used to compute waveforms for extreme 
mass ratio inspiral [3^, [4(| ■ 

To simplify the calculation we will use a normalized 
energy E r = E r /E, and define the following normalized 
impact parameter: 



< b < b* 



b + b c - 2b* 



b* <b<b r 



(32) 
(33) 



Though at energies probed by the LHC, and given that current 
experiments suggest the Planck scale cannot be very far below 
maximum LHC energies, the effects of spin and charge would be 
important to consider |38H . 



b c is a cut-off value for the impact parameter beyond 
which ([9|) ceases to offer a good approximation to the 
scattering behavior; geodesic integrations suggest b c w 
2b* in most cases. With the above normalization < 
6 < 2, and ([9]) becomes 



n(b) = -7 m |1 - b\. 



(34) 



In the above we have added the additional approximation 
that 77.(6 = 2) = 0; strictly speaking n is only zero at 
6 = (for non-spinning black holes) and in the limit as 
6 — * 00. We assume that the energy E r emitted during 
the process is simply a function of 6, and hence n by the 
above relation 



dE r /db - 



dE r dn 
dn db ' 



(35) 



and that a constant fraction e of the remaining energy is 
emitted per orbit in the whirl phase: 



dE r /dn = e(l - E r ) 
Integrating (f35|) with (I34|36|) then gives 



(36) 



E r (b) = 1- (1-E r0 ) (l-6) £ 
E r (b) = 1 - (6- l) e7 , 1 < 6 < 2 



< 6 < 1 



(37) 



For boundary conditions we have assumed all the energy 
is radiated for 6 = 1, E r0 is the energy radiated for the 
head-on collision case, and E drops to zero at 6 = 2. 

For interest, in the top panel of FigQT]we show a plot 
of (|37p with parameters plugged in from the low-speed 
system discussed in Sec lIIIl Of course, in this case E r (l) 
can at most be « .29 from the area theorem, though for 
values of 6 such that E r (b) <~ .29 ([37| is still a very 
good approximation. 

Returning to the ultra-relativistic problem, we now es- 
timate e using the quadrupole formula. For a circular 
orbit composed of equal mass point particles separated 
by a distance f and orbiting with angular velocity u>, 
where the over-bars (~) denote that the quantities have 
been normalized to the remaining energy 1 — E r in the 
system, the quadrupole formula gives 



47T 



f 4 u 5 . 



(38) 



So what values of f, uj and 7 to use? In the ultra- 
relativistic geodesic scattering problem, the geodesies ap- 
proach the light-ring (unstable circular photon geodesic) 
of the background spacetime at threshold. It also seems 
natural to guess that in this limit in the full problem 
the final angular momentum of the black hole will ap- 
proach extremality. The initial angular momentum of 
the system with critical impact parameter is J = b*E/2 
(restoring units); thus the initial effective Kerr parameter 
of the orbit is a a = 6*/2. A plot of 6* versus the back- 
ground Kerr parameter a for geodesic motion is shown 
in Fig[T2] — note that in the limit a/M — ► 1, do — > 1. To 
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estimate how a Q changes during evolution, we again use 



quadrupole physics, which says dJ/dt 



1 dE/dt, and 



define the instantaneous a D (t) = J{t)/E{t). We then get 



d{a /E) dE 1 



1 



dn 



2 ± 

dn E \ Euj E 2 



(39) 



In the limit a/M — > 1, J/E 2 — ► 1 initially, and using 
the Boyer-Lindquist value (fT5|) for u> for a photon on the 
light-ring in the extremal case, we get Elo — > 1/2. In 
other words, d(ao/E)/dn = in this limit, so at least 
to within the quadrupole approximation assuming an ex- 
tremal Kerr background for the geodesic analog in the 
ultra relativistic limit is self consistent. 

Though now we have a bit of a dilemma — in the ex- 
tremal limit there are no unstable circular orbits of Kerr, 
and hence 7 — + 00! In a practical setting (such as the 
LHC) one of course will not be at the limit, though given 
how sensitive 7 is to a approaching the limit it will be 
difficult to justify any crude estimates such as that out- 
lined in the previous paragraph to decide which value of 
a to use to determine b* and 7 from geodesic motion. 
In the bottom panel of FigfTT] we therefore plot several 
possibilities for E r (b) with a few values of a close to 1. 
We have used the limiting value of e qua d — * 7r/40, which 
is less sensitive to a than 7 in the limit. 

If FigfTT] gives a decent description of the ultra rel- 
ativistic particle scattering problem, then even though 
the cross section for black hole production will to a good 
approximation be 7r&* 2 /4, there can still be significant en- 
ergy loss to gravitational waves for b up to almost twice 
b* (or possibly even more, since recall that E r (b = 2) = 
was only an approximate boundary condition we im- 
posed), b* w 2E in the limit, which implies an effective 
cross section for measurable energy loss of AttE 2 . It is 
interesting that this number is identical to the order of 
magnitude estimate made by assuming the cross section 
is equal to ttR 2 , where R s — 2E is the Schwarzschild ra- 
dius corresponding to the initial center-of-mass energy. 
We finally note that a lower limit on the impact param- 
eter resulting in black hole formation can be computed 
by searching for trapped surfaces at the moment the two 
shock waves representing the r — * 00 particles collide; in 
[HI trapped surfaces were found for b <« 1.68E. 



VI. SUMMARY AND CONCLUSIONS 

In this paper we have described numerical simulations 
of the merger of a class of equal mass, zero initial spin, 
non-circular binary black hole systems in general rela- 
tivity. For a one parameter (k) family of solutions in- 
terpolating between a deflection of the binaries without 
merger at one extreme of the parameter, and a prompt 
merger at the other extreme, we provided evidence that 
there is a notion of a threshold of immediate merger at 
k ss k* during with the binary enters a tight near-circular 
whirl configuration before either separating or merging. 
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FIG. 11: Estimated fraction of the total energy radiated as a 
function of the impact parameter for the black hole scattering 
problem (|37[1 . The top panel illustrates a rest- mass dominated 
case, using parameters from the simulation results presented 
in Sec lIIII here cosmic censorship would place an upper bound 
for E r (l) of m 0.29. The bottom panel illustrates several 
curves from the kinetic energy dominated regime (note the 
different vertical scales in the two panels). For a concrete 
example we have used the value of 0.15 for Era, which is 
about half that of the upper limit computed using the trapped 
surface method (for a review of various methods see [4l[ ) . The 
order-of-magnitude calculation described in the text suggests 
that in the ultra relativistic limit the final state will be an 
extremal Kerr black hole (with negligible mass relative to the 
initial kinetic energy of the system) . In this limit the geodesic 
analogue calculation becomes very sensitive to the value of 
the final spin parameter a, and so we show several curves 
for different values of a. In the limit a — ► 1, 7 — ► 00, and 
the estimate E r (b) in (|37[) approaches a unit step function 
0(5/2). 
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FIG. 12: The critical impact parameter b* for equatorial 
geodesies with initial Lorentz F = 10 on a Kerr background 
with spin parameter a and mass m. 



The number of orbits n spent in the whirl is exponen- 
tially sensitive to the initial conditions: e" oc \k — fc*| 7 , 
where 7 is approximately a constant. The area theorem 
together with measurements of the energy lost to gravi- 
tational waves suggests this whirl behavior could persist 
for as many as 20-30 orbits for extremely fine tuned ini- 
tial data, though we have only been able to tune to w 5 
orbits due to limited computational resources. 

A second result of this work has been to show that 
similar behavior is observed in the analogous problem 
of the scattering of a one parameter family of equato- 
rial geodesies off a Kerr black hole. At threshold, the 
geodesic approaches one of the unstable circular orbits of 
Kerr at a radius r$ that depends on the particular fam- 
ily of geodesies. One quantitative similarity between the 
test particle and equal mass cases is we notice roughly the 
same scaling exponent 7 with geodesies orbiting a black 
hole with spin parameter a close to that of the black hole 
that forms in the merger case in the full problem, and ap- 
proaching an unstable orbit with radius similar to the co- 
ordinate separation of the black holes during their whirl 
phase. Comparison with the quadrupole formula for the 
emitted gravitational waves suggests the simulation co- 
ordinates are well adapted to the underlying physics. 

The close analogy between threshold geodesic scatter- 
ing and the one class of equal mass interactions studied 
here, together with an argument that such a threshold 
should exist for generic one parameter families of initial 
configurations with appropriate limiting cases, gives us 
some confidence in extrapolating this behavior to the ki- 
netic energy dominated regime. This is of relevance to 
parton scattering in the LHC if large extra dimension 



scenarios describe these interactions in the TeV regime, 
for then high energy partons could form black holes, and 
the scattering process would be well approximated by 
the collision of two black holes. We presented a sketch 
of how the geodesic analogue could be used to estimate 
the black hole formation cross section, and energy radi- 
ated to gravitational waves as a function of the impact 
parameter. At threshold it is conceivable that essentially 
all the kinetic energy is radiated as gravitational waves. 
Away from threshold significant amounts of energy could 
still be radiated even if a black hole does not form. 

Clearly, much of the above is pure speculation, though 
we believe is sufficiently interesting and relevant that it 
will be a fruitful endeavor to try to further establish (or 
disprove) the correspondence between geodesic scattering 
and the full problem. Simulating the ultra relativistic 
collision of black holes will computationally be a very 
challenging problem, and probably not possible to do 
beyond the 4-dimensional case without peta-fiop scale 
computing. Though if the correspondence could be es- 
tablish more strongly in the 4-dimensional case (which 
could be tackled with tera-fiop resources, at least for 
moderately relativistic energies), then comparable mass 
head-on collision simulations in higher dimensions to- 
gether with studies of geodesic scattering and estimates 
of energy and angular momentum emission in gravita- 
tional waves, similar to that performed in [44| for exam- 
ple, could be used to provide useful information for the 
higher dimensional case. 
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APPENDIX A: GEODESIC INTEGRATION 

In this appendix we describe the manner in which 
we integrate geodesies in the Kerr geometry. The orig- 
inal purpose of this geodesic integration method was 
for straight-forward incorporation into the GH evolution 
code, to study geodesic propagation in binary merger 
spacetimes. Thus the method does not take advantage of 
any of the symmetries of the underlying spacetime, nor 
uses advanced high-order ordinary differential equation 
(ODE) integrators. Nevertheless given the speed of con- 
temporary desktop PC's and that geodesic integration is 
a one dimensional evolution there is no problem in ob- 
taining sufficient accuracy for the purposes of the studies 
presented in the main sections of the paper. 



14 



Consider a spacetime with metric g a p 

ds 2 = g a pdx a dx^ , 



(Al) 



and a geodesic of the spacetime described via the para- 
metric representation of its curve x a = x a (X), where A 
is the affine parameter along the curve. The geodesic's 
tangent vector is 



a _ dx a (\) 
1 dX ~ 



(A2) 



defining prime (/) to denote differentiation with respect 
to A. For x a to represent a geodesic, u a must satisfy the 
geodesic equation: 



u' a + V = o, 



(A3) 



where T" s are the Christoffel symbols ([3]). In terms of 
the coordinate position along the curve x a (X), (|A3[) can 
be written as: 



r^x^x' 5 = 0. 



(A4) 



The causal character of the geodesic is given by the nor- 
malization of its tangent vector, as follows: 



u a u^g al3 



- 1 timelike 



u a u p g a fj = null 
u a vPg a p = 1 spacelike 



(A5) 
(A6) 
(A7) 



We are most interested in timelike and null geodesies at 
the moment. For timelike geodesies, the above normal- 
ization is equivalent to demanding that A measures the 
proper time of an observer moving along the geodesic. 
A does not have such a straight-forward interpretation 
for a null geodesic; in fact, we can re-parameterize any 
null geodesic via a linear transformation of the form 
s = aX + b, where a and b are constants. 



The geodesic equations (|A4I) are a set of four ODEs for 
the coordinate position of the corresponding "particle" as 
a function of affine time A. There are at least a couple 
of ways to proceed to solve these equations. One is to 
integrate (|A4|) directly as a set of second order ODEs, 
the other is to reduce it to a system of first order ODEs. 
We will take the former approach as it will fit "naturally" 
within the generalized harmonic evolution code metrics 
that we want to explore the geodesic structure of. An- 
other practical consideration in the code is that we need 
to integrate with respect to coordinate time t and not 
affine time A. Define x° = t, and let x k , k <G 1,2,3 be 
the three spatial coordinates, i.e. (a: 1 , x 2 , x 3 ) = (x,y, z). 
We thus want to solve for x k {t) for each geodesic. With 
the overdot (') denoting differentiation with respect to t, 
using the chain rule we get: 



x' k = xH 1 

c »* = X k t 12 + X k t" , 



Substituting (IA8IA9[) into !(A4[) gives 



x a t' 2 + x a t" + r^iVi' 2 = 0. 



Equation (|A4|) for t reads: 



t" + T^x' s = 0, 



(A10) 



(All) 



and using this, with (|A8[) again, the geodesic equation 
(|A10p for the spatial coordinates x k can be written as: 



x k + (T k s - x k T° jS ) x*x 5 = 0. 



(A12) 



This is the set of equations we will solve for timelike and 
null geodesies in a general spacetime g a p. 



1. Initial Conditions 

Since we are solving three, second order in time ODEs 
for each geodesic, we need six initial conditions per curve: 
the initial position x k (t — 0) and coordinate velocity 
x k (t = 0) of each particle. Note that we cannot choose ar- 
bitrary velocities, as our choices must be consistent with 
the normalization conditions (|A5IA6[) . There are several 
conceivable ways to ensure this; we will take the following 
route. 

We begin by choosing some arbitrary direction fc 7 = 
(0,k x ,k v ,k z ) for the curve to point in. We need the unit 
spatial vector in this direction; call it k 1 : 



Wtfg^s = 1 
and k 1 can be found by defining 

fc 7 = TVfc 7 



(A13) 



(A14) 



plugging this into the normalization condition (|A13|) , and 
solving for N . The tangent vector u a corresponding to a 
photon moving in the direction k 1 is then: 



u a — n a + k (null case) 



(A15) 



For a timelike curve, in addition to the direction k a we 
can choose a velocity v (< 1). The corresponding tangent 
vector in this case is 

u a = T (n a + vk a \ {timelike case), (A16) 

where T is the Lorentz gamma factor (here relative to an 
observer sitting at constant coordinate location and thus 
moving in the n a direction): 

r = -TpUr (A17) 



(A8) 

Now that we know u a , using (|A2[) and (|A8|) we can find 
(A9) the initial coordinate velocities of our geodesic curves. 
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2. Numerical Technique 

For compatibility with the metric evolution code we 
use a three time level scheme in a Cartesian coordinate 
system. The discrete version of the curve (x(t), y(t),z(t)) 
is represented by x l , y l , z l , where the time t l = iAt, with 
At the discretization scale. Time derivatives are com- 
puted using standard second order accurate stencils: 



/(*)!*=* 



i+1 



r 



f(t)\ t =t 



2At 

p+i _ 2 fi + f 
AP 



- 0{At 2 ) 
0{At 2 ). (A18) 



Written in terms of the residual Rj = Cj{q[) and the 
Jacobian matrix Jji = the resultant linear system 
that is solved for the correction S l q is: 



Jji ■ $q 



-Rj. 



(A22) 



After solving for the correction, the above steps are re- 
peated with the corrected solution serving as the new 
guess, and the iteration proceeds until the norm of the 
residual is below some specified tolerance. 

a. Second order accurate initial conditions 



Thus, when evolving from t = t l to t = t l+ , we require 
the locations of the geodesies at t — t % and t = t 1 ^ 1 . In 
the geodesic equation (|A12[) the Christoffel symbols are 
supplied as known functions of the spacetimc. 

When the time differences (|A18P are substituted into 
the geodesic equation (|A12[) for the time derivatives of 
the three coordinate positions of a geodesic, we end up 
with an algebraic system of three equations for three 
unknowns: (x l+1 , y l+1 , z l+1 ) . These equations are non- 
linear, and so an efficient method to solve for them is via 
Newton iteration. We write the system of equations as 



Cj{q k 



0, 



(A19) 



where j — 1,2,3 labels one of the equations, and we 
now use the notation q k , k = 1,2,3 to label one of the 
unknowns (i.e. q 1 — x l+1 ,q 2 — y l+1 ,q 3 = z l+1 ). Newton 
iteration then proceeds by writing the unknowns as a 
guess q k plus a correction 5q k , 



q k = q\ + Sq k 



(A20) 



linearizing about the guess, and solving for the correc- 
tions to first order: 

C 3 {q k ) = C,(q k + Sq k ) 

« Cjl^ + — f \ qk=4} ■ Sq l = 0. (A21) 



To begin evolving the geodesies at t = with a three 



time level scheme, we need the initial positions x , y 
at t = and the positions x~ \y~ 



,,0 ^0 



V- 1 



at t 



-At. 



One can use taylor expansions, the freely specifiable ini- 
tial conditions discussed in Sec. IA 1[ and the geodesic 
equations (|A12jl to initialize the past time values to the 
necessary level of accuracy. Use the subscript to refer 
to the analytic initial condition for a given variable; for 
example, for x: 



xo = x(t)\ t=0 , 
x = x(t)\ t=0 , 
x = x(t)\ t =o- 

Then, to second order 

x^ 1 — Xq — XoAt + Xq 



At 2 



(A23) 
(A24) 
(A25) 



(A26) 



The initial position x$ is freely specifiable, the initial 
velocity ±q is calculated as described in (|A 1[) . and the 
geodesic equations (|A12I) are used to solve for £§: 



(A27) 



Note the summation over all initial velocities in the 
above, including over t = 1. 
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